#!/usr/bin/perl -w
# use strict;
# use warnings;
################################################################
# Getting input parameters of preprocess.pl program      #
################################################################

my ($USAGE) =  "\n\n\t****** USAGE of preprocess.pl PROGRAM ******\n\n\n\n\tUSAGE: perl $0 <panam.ini file> \n\n\n\n";
die "$USAGE" if ((scalar(@ARGV))< 1);
my $option_file = $ARGV[0];
chomp($option_file);

die "\n\n\t=> Cannot find configuration file: $option_file.\n\n" unless (-e $option_file);
die "\n\n\t=> Configuration file, $option_file, appears to be empty!\n\n" if (-z $option_file);
open(OFILE, "<$option_file") || die "Cannot open input file : $option_file\n";

my $inputSeqName;
my $inputQualName;
my $lengthCut;
my $qualityCut;
my $NGS_id_Results;
my $match;
my $forward;
my $reverse;
my $pmismatch ;
my $Nbase; #### ajout 31/5

#my $fuzznuc=`which fuzznuc`; 
#chomp ($fuzznuc);
my $fuzznuc = "/home/panam/panam_v3/bin/EMBOSS-6.5.7/emboss/./fuzznuc" ;

my $path;

while (<OFILE>) {
	chomp;
	my $line = $_;

# discard blank line
	if ($line =~ m/^\s*$/g) {
		next;
	}

# discard comment line
	elsif($line =~ m/^\s*#/g) {
		next;
	}

# get input options line
        else {
	
		if($line=~m/INPUT_FILE/gi) {
			$line=~s/INPUT_FILE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$inputSeqName = $line; #print "$inputSeqName\n";
			}
			chomp ($inputSeqName);
			# See if the 454 Sequence file exist and contains something
			die "\n\n\tCannot find 454 sequence file: $inputSeqName. Check panam.ini.\n\n" unless (-e $inputSeqName);
			die "\n\n\t454 sequence file, $inputSeqName, appears to be empty!\n\n" if (-z $inputSeqName);

			my $ligne1= `sed -n '1p' $inputSeqName`;
			my $ligne2= `sed -n '2p' $inputSeqName`;
			##### modif 15/5/2012
			if (($ligne1 !~ /^>/) or ($ligne2 !~ /^[ATCGUatcgu]/)) {
				die "\n\n\tYour file $inputSeqName does not seem to contain sequences in fasta format.\n\n";
			}
		}

		if($line=~m/QUALITY_FILE/gi) {
			$line=~s/QUALITY_FILE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$inputQualName = $line; #print "$inputQualName\n";
			}
			chomp ($inputQualName);
			# See if the 454 Quality file exist and contains something
			die "\n\n\tCannot find quality file: $inputQualName. Check panam.ini.\n\n" unless (-e $inputQualName);
			die "\n\n\tQuality file, $inputQualName, appears to be empty!\n\n" if (-z $inputQualName);

			my $ligne1= `sed -n '1p' $inputQualName`;
			my $ligne2= `sed -n '2p' $inputQualName`;	
			if (($ligne1 !~ /^>/) or ($ligne2 !~ /\d/)) {
				die "\n\n\tYour file $inputQualName does not seem to be in the right format.\n\n";
			}
		}

		if($line=~m/LOWER_SEQUENCE_LENGTH_CUTOF/gi) {
			$line=~s/LOWER_SEQUENCE_LENGTH_CUTOF//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$lengthCut = $line; #print "$lengthCut\n";
			}
			die "\n\n\t Length cutoff missed. Check panam.ini.\n\n" unless ($lengthCut ne "");
		}

		if($line=~m/AVERAGE_CUTOFF_QUALITY_VALUE/gi) {
			$line=~s/AVERAGE_CUTOFF_QUALITY_VALUE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$qualityCut = $line; #print "$qualityCut\n";
			}
	 		die "\n\n\t Quality cutoff missed. Check panam.ini.\n\n" unless ($qualityCut ne "");
		}
	
		if($line=~m/PMISMATCH/gi) {
			$line=~s/PMISMATCH//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$pmismatch = $line; #print "$lengthCut\n";
			}
			die "\n\n\t PMISMATCH should be a percentage value! Check panam.ini.\n\n" unless (($pmismatch >=0 and ($pmismatch <= 100)));
		}

		if (!(defined $pmismatch)) {$pmismatch = 0}
		
		if($line=~m/454_RUN_IDENTIFIER/gi) { ###### dossier !
			$line=~s/454_RUN_IDENTIFIER//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$NGS_id_Results = $line; 
			}
	 		die "\nOutput directory for NGS analyses is not defined. Check panam.ini.\n\n" unless ( defined $NGS_id_Results);
			unless (-d $NGS_id_Results) { mkdir $NGS_id_Results || die "Could not create Directory $NGS_id_Results $!\n"; }
		}

		if($line=~m/MATCH_PRIMER/gi) {
			$line=~s/MATCH_PRIMER//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$match = $line;  
			}
	 		die "\n\n\t Match primer value missed. Check panam.ini.\n\n" unless ($match ne "");
			die "\n\n\t Match primer value is not correct. Check panam.ini.\n\n" unless (($match eq "both") or ($match eq "forward") or ($match eq "reverse"));
		}

		if ($line=~m/FUZZNUC_PATH/gi) {
			$line=~s/FUZZNUC_PATH//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi; 
			if($line) {
				$path = $line; 
			}
		}
		
		if($line=~m/FORWARD_PRIMER_SEQUENCE/gi) {
			$line=~s/FORWARD_PRIMER_SEQUENCE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$forward = $line; 
			}
		}
		
		if($line=~m/REVERSE_PRIMER_SEQUENCE/gi) {
			$line=~s/REVERSE_PRIMER_SEQUENCE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$reverse = $line; 
			}
		}

		###########ajout 31/5
		if($line=~m/AMBIGUOUS/gi) {
			$line=~s/AMBIGUOUS//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi; 
			if($line) {
				$Nbase = $line; 
			}
			die "\n\n\t Ambiguous base value missed. Check panam.ini.\n\n" unless ($Nbase ne "");
			die "\n\n\t Ambiguous base value should be yes or no. Check panam.ini.\n\n" unless (($Nbase eq "yes") or ($Nbase eq "no")) ;
	  
		}
		########### fin ajout 31/5
	}
}

die "\n\n\t Match primer value is set at 'both' but can not find forward primer. Check panam.ini.\n\n" if (($match eq "both") and !(defined $forward));
die "\n\n\t Match primer value is set at 'forward' but can not find forward primer. Check panam.ini.\n\n" if (($match eq "forward") and !(defined $forward));
die "\n\n\t Forward primer value is not correct. Check panam.ini.\n\n" if ((defined $forward) and ($forward!~ /[ATCGUatcgu]/));

die "\n\n\t Match primer value is set at 'both' but can not find reverse primer. Check panam.ini.\n\n" if (($match eq "both") and !(defined $reverse));
die "\n\n\t Match primer value is set at 'reverse' but can not find reverse primer. Check panam.ini.\n\n" if (($match eq "reverse") and !(defined $reverse));
die "\n\n\t Reverse primer value is not correct. Check panam.ini.\n\n" if ((defined $reverse) and ($reverse!~ /[ATCGUatcgu]/));


if (!(-e $fuzznuc)) {
	if (defined $path) {
		$fuzznuc=$path."/fuzznuc" ;
	}
}
if(!(-e $fuzznuc)) {
	die "\n\n\tfuzznuc cannot be found\n";
}

############################################ trimming sequences ######################################################""

print "\n\tTrimming sequences ...\n\n"; 

if (-d $NGS_id_Results."/tmp") {
	`rm -r $NGS_id_Results/tmp`;
	`mkdir $NGS_id_Results/tmp`;
}
else {`mkdir $NGS_id_Results/tmp`}

my $outputName;
$outputName = "tmp/sample_all_oldIds_tmp";
open(INPUTSEQ, $inputSeqName);
open(INPUTQUAL, $inputQualName);
open OUTPUT, ">$NGS_id_Results/$outputName" or die $!; 

my $Rejected = -1;
my @FinalTrim = ();
my $end = 0;
my $max = 0;
my $sum = 0;
my $start = 0;
my $first = 0;
my $lineNum = 0;


while(my $lineSeq = <INPUTSEQ>)
{
	my $lineQual = <INPUTQUAL>; $lineQual=~s/^\s+//;
	if($lineSeq =~ />/)
	{
		my $length = ($end + 1) - $start;
		if($length < $lengthCut)
		{
			$Rejected = 1;
		}
		if($Rejected == 0)
		{
			print OUTPUT "$header"; 
			print OUTPUT "\n";
			$countDown = 60;
			for($a = $start; $a <= $end; $a++)
			{
				$countDown--;
				if (defined $FinalTrim[$a]) {print OUTPUT "$FinalTrim[$a]";}
				if($countDown == 0)
				{
					print OUTPUT "\n";
					$countDown = 60;
				}
			}
			print OUTPUT "\n";
		}
		$Rejected = 0;
	
		#Trim off heading

#############################
		if ($lineSeq=~ / /) {$space = index($lineSeq, " ") + 1}
		elsif ($lineSeq=~ /\n/) {$space = index($lineSeq, "\n") + 1}
#############################	  

		$header = substr($lineSeq, 0, $space);
		chomp($header);
		@FinalTrim = ();
		$max = 0;
		$sum = 0;
		$first = 0;
		$lineNum = 0;
	}
	elsif($Rejected == 0)
	{
		#Store for use when trimming bases
		@bases = split(//, $lineSeq);
		$size  = scalar @bases;
		$size--;
		for($a = 0; $a < $size; $a++)
		{
			push(@FinalTrim, $bases[$a]);
		}
	
		chomp($lineQual);
		@line = split(/ /, $lineQual);
		$size = scalar @line;
		for($a = 0; $a < $size; $a++)
		{ 
			$sum += ($line[$a] - $qualityCut);
			if($sum > $max)
			{
				$max = $sum;
				$end = $a + (60*$lineNum);
				$start = $first;
			}
			if($sum < 0)
			{
				$sum = 0;
				$first = $a + (60*$lineNum);
			}
		}
		$lineNum++;
	}
}

close INPUTSEQ;
close INPUTQUAL;
close OUTPUT;

open (R, $NGS_id_Results."/tmp/sample_all_oldIds_tmp") || die "can not open file";
my %seq; my $a; 
my %corresp; my $ii=1;
my $id;

############## ajout 28/5
my %ident;
###############fin ajout 28/5

while(<R>){
	if ($_=~ /^>(.*?)\s+/) {
		$a =$1; 

		##########ajout 28/5
		if (exists $ident{$a}) {
			 $a.="_".$ii;
		}
		else {
			$ident{$a} = 1
		}
		######### fin ajout 28/5

		###########ajout1 pour correspondance
		$id = "seq".$ii; 
		$corresp{$id} = $a;
		$ii++;
		#######"fin ajout1
		
	}
	else {
		$seq{$id}.=$_
	}
}

##############ajout2 pour correspondance
# open (T, ">".$NGS_id_Results."/seq_correspondance");		## VERSION 3.01 ## inutile
open (RR, ">".$NGS_id_Results."/tmp/sample_all_tmp") || die "can not open file";
foreach my $mp (keys %corresp) {
#	print T "$mp\t$corresp{$mp}\n";				## VERSION 3.01 ## inutile
	print RR ">$mp\n$seq{$mp}\n"
}
# close T;				## VERSION 3.01 ## inutile
close RR;
##############fin ajout2

`$fuzznuc -sequence $NGS_id_Results/tmp/sample_all_tmp -pattern $forward -complement Y -pmismatch $pmismatch -outfile $NGS_id_Results/tmp/outputF_tmp.fuzznuc -rformat excel` ;
`$fuzznuc -sequence $NGS_id_Results/tmp/sample_all_tmp -pattern $reverse -complement Y -pmismatch $pmismatch -outfile $NGS_id_Results/tmp/outputR_tmp.fuzznuc -rformat excel` ;

my %fuzz_F; my %fuzz_R;

open (FUZF, $NGS_id_Results."/tmp/outputF_tmp.fuzznuc") ;
open (RFF, ">".$NGS_id_Results."/tmp/outputF_Ids_tmp.fuzznuc") ;
while (<FUZF>) {
	if ($_=~ /SeqName/) { print RFF "$_" }
	else {
		my @t = split (/\t/, $_); 
		$fuzz_F{$t[0]}++;
		print RFF "$corresp{$t[0]}\t$_";
	}
}	
close FUZF;
close RFF;

open (FUZR, $NGS_id_Results."/tmp/outputR_tmp.fuzznuc") ;
open (RFR, ">".$NGS_id_Results."/tmp/outputR_Ids_tmp.fuzznuc") ;
my %h;
while (<FUZR>) {
	if ($_=~ /SeqName/) { print RFR "$_"}
	else {
		my @t = split (/\t/, $_);
		$fuzz_R{$t[0]}++;
		print RFR "$corresp{$t[0]}\t$_";
	}
}	
close FUZR;
close RFR;


my $i ;
if ($match eq "forward") {
	$i=0;
	open (FF, ">".$NGS_id_Results."/sample_all") || die "can not open file";
	foreach my $k (keys %seq) {
		if (defined $fuzz_F{$k}) {
			if ($Nbase eq "yes") {
				if ($seq{$k} !~ /N/i) {
					$i++;
					print FF ">$corresp{$k}\n$seq{$k}";
				}
			}
			if ($Nbase eq "no") {
				$i++;
				print FF ">$corresp{$k}\n$seq{$k}";
			}
		}
	}
}


elsif ($match eq "reverse") {
	$i=0;
	open (FF, ">".$NGS_id_Results."/sample_all") || die "can not open file";
	foreach my $k (keys %seq) {
		if (defined $fuzz_R{$k}) {
			if ($Nbase eq "yes") {
				if ($seq{$k} !~ /N/i) {
					$i++;
					print FF ">$corresp{$k}\n$seq{$k}";
				}
			}
			if ($Nbase eq "no") {
				$i++;
				print FF ">$corresp{$k}\n$seq{$k}";
			}
		}
	}

}


elsif ($match eq "both") {
	open (FF, ">".$NGS_id_Results."/sample_all") || die "can not open file";
	foreach my $k (keys %seq) { 
		if ((exists $fuzz_F{$k}) and (exists $fuzz_R{$k})) { 
			if ($Nbase eq "yes") {
				if ($seq{$k} !~ /N/i) {
					$i++;
					print FF ">$corresp{$k}\n$seq{$k}";
				}
			}
			if ($Nbase eq "no") {
				$i++;
				print FF ">$corresp{$k}\n$seq{$k}";
			}
		}
	}
}

#####################"" ajout 28/5 (if $match eq "none")
else {
	open (FF, ">".$NGS_id_Results."/sample_all") || die "can not open file";
	foreach my $k (keys %seq) { 
		if ($Nbase eq "yes") {
			if ($seq{$k} !~ /N/i) {
				$i++;
				print FF ">$corresp{$k}\n$seq{$k}";
			}
		}
		if ($Nbase eq "no") {
			$i++;
			print FF ">$corresp{$k}\n$seq{$k}";
		}
	}
}
##########################"" fin ajout 28/5

# 
open (L ,">$NGS_id_Results/PREPROCESS_LOG");
print L"Processed file: $inputSeqName\nQuality cut: $qualityCut\nLength cut: $lengthCut\nForward primer: $forward\nReverse primer: $reverse\nMismatch value: $pmismatch\n";

print "Successful trimming, You have $i retained sequences, they are in the file sample_all.\nYou can use this file for further preprocessing.\n\n";

#`rm -r $NGS_id_Results/tmp`;
